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ABSTRACT 


It is shown that the Flux-Corrected Transport (FCT) 
algorithm correctly computes the propagation velocity and 
Chapman-Jouguet pressure for a laser-supported detonation 
(LSD) wave, even for zone size large compared to the 
radiation absorpticn length. A two-dimensional FCT computer 
program for analysis of LSD wave propagation phenomena is 
presented. A trajectory method for analysis of two-dimen- 
sional laser beam propagation in the paraxial Fresnel ap- 
proximation, with continuously varying refractive index, is 
described. 


1. INTRODUCTION 


The developments described below constitute the 
initial portion of an effort to develop and apply efficient, 
versatile analytic techniques for the analysis of laser beam 
propagation and laser-target interactions. The two problems 
addressed in the report are described below. 


1.1 EVALUATION OF NON-ADIABATIC SHOCK PROPAGATION 


A two-dimensional hydrodynamics computer code was 
written using the FCT (Flux-Corrected Transport) method de- 
veloped at NRL [References 1-5]. The performance of this 
code in describing the propagation of laser-supported 
detonation waves was studied in some detail from the stand- 
point of thermodynamic validity for the case of one-dimen- 
sional wave propagation; initial testing of two-dimensional 
detonation wave propagation was also carried out. The re- 
sults in both cases appear to be highly satisfactory. The 
code listing is presented in Appendix I. 


1.2 LASER-BEAM PROPAGATION IN PLASMAS 


Two-dimensional methods for computing refraction and 
Fresnel diffraction of laser beams in plasmas and partially- 
ionized gases were devised, using a trajectory approach of 


the general type suggested by Glass. !! A computer code em- 


ploying these methods is under development. A code listing 
is presented in Appendix II, representing the implementation 
of the text formulation at an incomplete stage of testing. 
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2. FLUX-CORRECTED TRANSPORT (FCT) APPLIED TO ONE-DIMENSIONAL 
LASER-SUPPORTED DETONATION (LSD) WAVE PROPAGATION 


In this section the results of a series of calcula- 
tions using the FCT method to describe the propagation of one- 
dimensional LSD waves are presented. The object is to deter- 
mine the accuracy with which the numerical method reproduces 
such details of a steady-state wave as the spike pressure, 
the Chapman-Jouguet (CJ) pressure, and the propagation i 
velocity. | 


Zot EQUATION OF STATE 


The equation of state used for these calculations is 
that of an ideal gas modified to allow for single ionization. 
The specific internal energy is related to the temperature 
and specific volume by 


ge» 2 ee E + 3x + a | (1) 
2 kT : 
where R is the gas constant (per gram), T is the temperature 
in kelvins, f is the number of degrees of freedom, k is the 
Boltzmann constant, Ey is the ionization potential, and x is 
the fraction of molecules which are ionized. It is related 
to the temperature and specific volume, V, by the Saha 


equation 


x ae ayp3/2 


I-x exp (- €./kT) (2) 


where the constant A is given by 


5.93 M_ =|" 


h 


’ 
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where M is the gram molecular weight, No is Avogadro's num- 
ber, m, is the electronic mass, and h is Planck's constant. 
With 3. = 2 and 3;/9, = 1.9 (the weight factor for ions 
relative to atoms), the numerical value for A corresponding 
to a molecular weight of M = 28 g/mole is 

A = 4.2662 x 10°’ gem x °/?, 

Given values for the specific volume and internal 
energy (the calling arguments from the hydrodynamic sub- 
routine), the temperature is obtained from Equations (1) and 
(2) by an iterative procedure. The pressure is then calculated 
from the ideal gas relation, 


PV = R(1l + x)T. (3) 


The equation of state subroutine used for the above 
calculations also returns values of the sound speed for use 
in determining the time-step. A Fortran listing is given in 
Appendix I. 


Zee STRUCTURE OF A ONE-DIMENSIONAL LSD WAVE 


A qualitative sketch of the structure of a steady- 
state, one-dimensional LSD wave is given in Figure l. An 
initial shock at the spike pressure preheats and partially 
ionizes the gas, providing the free electrons necessary to 
initiate absorption from the incident beam. As the energy 
absorbed increases through the absorption front, the pres- 
sure decreases and the temperature and specific volume in- 
crease until the CJ point is reached, at which point 
essentially all of the incident beam energy has been absorbed 
by the gas. 


The quantitative determination of the spike and CJ 
pressures for a given laser flux is indicated in Figure 2. 
The Hugoniot relation for conservation of energy for a steady 
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Figure 2. Hugoniots for an ionizable ideal gas corresponding 
to zero energy absorption and to complete energy 
absorption for a laser beam flux of 108 watt/cm2. 


state absorption wave connecting an initial state Po? ie E 


fo) 
to a final state P, V, E (Figure 3) is 


me! ae G 
E Ey = 5 (P + Fo) Vv) + BoD r (4) 
where G is the laser flux (erg/cm2/sec) and D is the wave 
velocity, which for a steady~state wave is related to the 


pressure and volume by 


P- Po 

Rees °. (5) 
fe) Le V 

The curve labeled "Complete Absorption" in Figure 2 

is the Hugoniot as calculated from Equation (4), the P, V, E 
equation of state as described in Section 2.1, and a laser 
flux G of 108 watt/em". The initial state is Pe = 10° dyne/ 
cm’, T, = 298.16 kelvins, <a e 885.4 em?/g,and E, = 2.213 x 


10? erg/g. The parameters used in the equation of state are: 


f =% 
R = 2.9694 x 10° erg/g/K 
Ey = 14.4 ev 


A = 4.2662 x 10°’ g om? «77/2 

The CJ point is that point at which the Rayleigh line 
is tangent to the Hugoniot curve, as indicated in Figure 2. 
The spike pressure is given by the intersection of this same 
Rayleigh line with the Hugoniot corresponding to zero energy 
absorption (calculated from Equation (4) with G = 0). For 
this equation of state and a flux of 108 watt/cm“, the spike 
pressure is 1.0 kbar, the CJ pressure is 0.45 kbar, and the 
wave speed is 10.0 * 10° cm/sec. 


eeg§——_———__— Flux G (erg/em*/sec) 


-_s 
ie) 


(Poot e) 


Absorption Wave 


Figure 3. Steady-state absorption wave connecting initial 
state P ae Vint Eo to find state P, V, E. 
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PAS NUMERICAL METHOD 


The Fortran listings of the one-dimensional FCT modules 


as given by Boris, 1] 


modified where necessary to allow for 

| running on the CDC-7600 and Univac 1108 computers, have been 
incorporated -nto a preliminary two-dimensional test code. 

A complete listing is given in Appendix I. The subroutine 
UPDATE is a realization of the two-dimensional time-splitting 
procedure as outlined in Reference [2]. The calculations 
reported here were done using the one-dimensional option, 
which is capable of operating in either the Eulerian or the 
Lagrangian mode. 


The flux-corrected transport method has been discussed 


[I=5] 


in a series of papers. Briefly, it solves the general- 


ized continuity equation of the form 


(x™ gu) = Source Terms , (6) 


where the terms on the left are the one-dimensional equiva- 
lent of (36/3t) + V * (¢u). Specifically, in the present 
case the algorithm is applied in succession to the three 


equations, 
—~ a (pu) = 0, (7) 
2 (ou) + & (pu?) = - , (8) | 
G64 2 (eu) = - & (Pu), (9) 
where 
ey, (= + 5 u’ | . (10) 
9 


R-3564 


An alternative form of the last equation is 
S— (cE) + 2= (cka) = = p 22 (11) 
dt ox 3 


In the above equations, which represent the conservation of 
mass, momentum and energy, respectively, 0 is the density, u 

is the material velocity, P is the pressure, x is the Eulerian 
space coordinate, and E is the specific internal energy (erg/g). 


2.3.1 Absorption Coefficient 


The laser beam intensity as a function of position for 
one-dimensional propagation in an absorbing medium is given 
by 


= expe 


x 
i k dx 
x 


fe) 


The resulting expression used to compute the increase in 
specific internal energy within a computational zone in time 
St is 


‘ e bs At 3 
Oy SE, * Tay F exp ( K,, jax l)] 38 erg/cm’, (12) 


where I is the intensity (erg/cm2/sec) incident on the 


m+1 
right interface of the zone. The absorption coefficient k,, 
for a laser frequency v is given by the sum of the absorp- 
tion coefficients for free-free transitions in the ionic 


and atomic fields, 
K 2k. ;, *£& ’ (13) 
[7] 


which are given by the formulae 


10 


=172 


(ore fl = exp(- busy) 82.4) 


Ba 9 3 x (l-x) 13/2 (7m?) {1 - exp(- hv/kT)] (15) 


2 


where Z is the degree of ionization (equal to 1.0 for these 
calculations), X is the laser wavelength in um, x is the 
fraction of ionization from Equation (2), M is the gram 
molecular weight, and o is the collisional cross-section be- 
tween electrons and molecules, which was assigned the value 

0 cm. The expressions for the numerical coefficients in 
the above expressions are 


6,_2 
ape ee a 
1 3m_k q 
e mc 
e 
= 4.97 x 1072 microns? Kl/2 em” ote 
2-2 
és (2 Se eee 
2 ™ he? 
« 4.69 « 10°! wdieren”™” K°°’* em” aote™* 


Actually, the precise form of the density and tempera- 
ture dependence of the absorption coefficient was of little 
consequence for these calculations, since the detailed 
structure within the absorption wave was not of interest. 
For example, at the CJ point corresponding to a flux of 10 
watt/om?, the specific volume is 524 cm?/gm, the temperature 
is 41,500 kelvins, and the absorption coefficient as cal- 


8 


culated from Equations (14) and (15) for a 10.6 micron 


laser wavelength is 4100 em -, i.e., the absorption front 


is of the order of i cm thick. Thus, extremely fine zoning 


rasan abate aa ececasae eae ee ae 
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is required to resolve the structure of the absorption front 
itself. 


2.3.2 Test Calculations 


The ability of the FCT method to describe the develop- 
ment and propagation of an LSD wave using a zoning which is 
coarse with respect to the physical thickness of the absorp- 
tion front is shown in Figure 4. The calculational con- 
figuration consists of 200 zones, each 0.01 cm thick, with 
a rigid wall for the left boundary. The wave was "ignited" 
by an initial heating of the first three computational 
zones. An initial energy in these zones of 3.96 x 40> 
erg/g at the ambient density of 1.1295 x to7° g/em* results 
in an initial pressure of 96.8 bars, and an initial tempera- 
ture of 20,360 kelvins, which produces sufficient initial 


absorption for wave development. The energy in this cal- 
culation was updated with Equation (9). 


Since the zone width of 0.01 cm is approximately 40 
times the absorption length at the CJ point, all of the beam 


energy is absorbed in one computational zone, and, as ex- 
pected, such details as the 1.0 kbar spike pressure are not 
reproduced by the calculation. The other features of the 
wave profile are reproduced quite well, however. As shown 
in Figure 4, the wave profiles behind the "pseudo" spike 
correlate very well with the analytical value of 0.45 kbar 
at the CJ point and with the analytical values of the slopes 
at the CJ point as obtained with a simple rarefaction wave 
centered at x = 0: 


aP Cc 

come = (16) 

(#) 2 

CJ p atc-) 
F + 5 eo x 


where the square of the sound speed, c*, and the density, 0p, 


12 
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are the values at the CJ point, 


co? = 3.56 x 107+ em*/sec* 
= 1.91 x 107? g/cm? 
2 
ais ) = 271 (cm*/sec”) / (dyne/cm?) i 


In the interval between 100 cycles (0.466 usec) and 400 cycles 
(1.668 usec), the peak wave pressure progresses from zone 47 
to zone 167. The corresponding average wave speed D is then 
9.98 x 10° 
10.0 x 10° cm/sec. Finally, the steady-state pressure at 
the left boundary is about 0.16 kbar, also in agreement with 


cm/sec, in agreement with the analytical value of 


the analytic value of v 0.159 kbar. The conclusion is that 
the FCT method can propagate an LSD wave properly even with 
zoning which is much too coarse to resolve the structure with- 
in the actual absorption front. 


The results of a calculation identical to that de- 
scribed above except that the zone size has been decreased 
from 0.01 cm to 5.0 x i cm are given in Figure 5. For 
this calculation the energy absorption front is about five 
zones thick, and as a result, the calculated spike pressure 
of v 0.9 kbar is more nearly equal to the analytic value of 
1.0 kbar. The ripples superimposed on the rarefaction fol- 
lowing the spike are not physical, but are an example of a 
nonlinear numerical error called "terracing". This type of 
error, which occurs on the flanks of steep gradients, is 
discussed in Reference [2]. 


The two calculations described above were done using 
the FCT modules in the Eulerian mode, i.e., the zone bound- 
aries were fixed in space. The results of a calculation 
using the Lagrangian mode are shown in Figure 6. The initial 


14 


R-3564 


*yzbueT uotydzosqe ey JO TepzIO 9Yz JO ST WO ¢_OT x 0°S JO YRPTM ouOZ 


* uo/33eM got 


au0zZ 
002 OST OoT OS 0 


15 


eansseig 
LO 


(zeqy) eanssezg 


SeTOAD 00S 


———— 


R-3564 


*yqbuetT 
uotzdiosqe ay} 03 3Oedser YIM BSTROD ST WO TO°Q JO YIPTM sUOZ TeTRTUT 


* 7uo/33eM gO JO XNTJ & TOF BAEM CST & FO uUOTReTNOTeO LOI uethuezbey °*g9 oanhty 


9uoZ 


002 OST OOoT 0S 0 


(zeqy) eanssezg 


seTOAD 00S 


aedots [erz#tuLl pue eansserg CD TROT Ateuy &* 


16 


R-3564 


zone width was 0.01 cm, as in the first calculation, but the 
zone boundaries were allowed to move with the material 
velocity at each point. This calculation also reproduced the 
Main features of the wave quite well, even though the initial 
zone width was very large compared with the absorption length. 
The average wave speed from cycle 100 (0.3999 wsec) to cycle 
500 (1.8182 usec) was 10.0 * 10° cm/sec, in agreement with 
the analytic value of 10.0 x 10° em/sec for a flux of 108 
watt/em-. The wave profiles also correlate well with the 
pressure and the initial slope at the CJ point. It will be 
noticed from Figure 6 that the peak pressure of the pseudo 
spike fluctuates considerably as the wave progresses. This 
fluctuation is a result of the coarse zoning, but it is not 
an instability; it arises because at times the conditions 

at the shock front are such that the energy absorption occurs 
across two computational zones instead of one. When this 
happens, the pressure in the first zone approaches more nearly 
the analytical spike pressure. 


The results of a calculation to check the use of an 
artificially long absorption length are shown in Figure 7. 
The zone width for this Eulerian calculation was 0.01 cm; but 
if the value of k, as calculated from Equations (14) and (15) 
was such that k [Ax] > 1.0, it was reset to 1.0 for that 
zone. In other words, the beam intensity was not allowed to 
decrease by more than a factor of e- across any computa- 
tional zone. The result of this numerical limiting is an 
absorption-front thickness of four or five computational 
zones; and, as expected, the peak spike pressure is close 
to the analytical value of 1.0 kbar for a flux of 108 watt/ 
cm?, A disadvantage of artificially limiting the absorption 
coefficient is that the time required to attain steady i 
state is unrealistically long. From Figure 7 it is seen that 
even after 400 computational cycles, or a travel of 1.4 cm, 
the pressure behind the spike has still not reached the steady- 
state CJ value of 0.45 kbar. 
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2.4 CONCLUSIONS | 


The FCT method is a convenient and satisfactory 
numerical procedure for describing the propagation of one- | 
dimensional LSD waves. In particular, the method yields 
satisfactory results even for zoning which is much too 
coarse to resolve the structure within the energy absorption 


| 
front. 
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3. LASER BEAM PROPAGATION IN GASES AND PLASMAS 


The propagation of a laser beam is a highly special- 
ized problem from the standpoint of radiation transport in 
general, since the radiation field is very sharply peaked 
with respect to both angular distribution and frequency dis- 
tribution. These characteristics can only be maintained in 
a material medium whose properties vary continuously and 
Slowly over finite intervals of space and time, so that the 
variations may be considered as perturbations within these 
intervals. At the boundaries of such intervals, discon- 
tinuous material properties may be taken into account in 
the beam model by introducing additional beams; such may 
also originate in the interior of an interval by means of 
stimulated Raman scattering and similar processes. In this 
section the propagation of a single beam within a continuous 
medium will be discussed. The properties of the medium will 
be assumed given, although they may depend in various ways 
upon the radiation field itself. 


For a beam with axis in the z direction in a medium 
with complex refractive index, n, the amplitude E satisfies 
the scalar wave equation 

2 2 «2 


e+ ib - 2 28-0 (17) 


az Cu cae 


where the V operator is restricted to the x-y plane. The 
perturbation formulation discussed above suggests that the 
amplitude be written as the product of a slowly-varying 
envelope and a fast unperturbed part: 


i(k,z - wt) 


E(x,¥,2,t) _ W(x,y,z,t)e ° (18) 


where Ws, is the central laser beam frequency and ky Now /cr 


20 
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with ny being a constant real index of refraction assumed as 
a base of reference. Let 


w | 
a= 2 = Imn = Imk (19) 
Re n - ny Re k - ky 
i ee (20) 
Mo Ko 


represent the local absorption coefficient and index pertur- 
bation, respectively, then on dropping terms of second order 
in d/dz, 3/dt, a, u, and oe there results the time-depen- 
dent paraxial wave equation 


3 


<< 


2 : 1 
Vou + 2ik, ( + = 


ale 


)+ (ix,0 rs 2x21) i=* (21) 


Q 
N 


[8] 


directly for the real and imaginary parts of y in the time- 
[6] followed the 
alternative approach of separating into amplitude and phase 


where v is the group velocity. Fleck solved this equation 


independent, axially symmetric case. Glass 


parts. The latter method leads, as will be shown, to a 
formulation in terms of variables which have simple geometri- 
cal interpretations and are well suited to calculations with 
adaptive coordinate grids. We have adopted this latter 
approach in the present effort. 


Following Reference [6] we introduce real intensity 
and inclination variables I and u by the relations 


y = Fe’? (22) 
u=> V6 (23) 
(eo) 


With these definitions, the real and imaginary parts of (21) 


yield the pair of equations 


 h—h—h—a————— ee 
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3 Lo 
($e +See tu v) I=-(0+V-u) I (24) 
ee ae : o 
(Geteaets Rp ee Pin eB) (25) 
where the diffraction potential P is defined as 

2 2 
oes Yor. £ 7¥s 
P= 5 (SE | (26) 


These equations may be integrated along characteristics, or 
trajectories, which are the analogs in Fresnel diffraction 
theory of the rays of geometric optics. These trajectories 
are defined by the pair of equations 


ru (27) 


ah 
€ = = (28) 
where r is the radius vector projection in the plane normal 
to the z axis and t is the time. Along such trajectories 


(24) and (25) take the form 


fe= (+ * uy I (29) 


V(p + P) (30) 


Te) 
i} 


Let us consider an arbitrary portion of the beam cross-sec- 
tion for given z, t. This is a region R of the xy plane, 
bounded by a set of points r at which (27) holds. On inte- 
grating (24) over this region, with boundary terms defined 
by (27), one finds 


22 
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i a7 d : 
Lf (a1 + = at) ax dy + ff I dx dy = 0 (33) 
R R 


which implies that energy flows only along, not across, tra- 
jectories. We may, in fact, replace (29) by (31) in order 
to have an explicit equation for energy conservation. Simi- 
larly, with the aid of (27) we may replace (30) by 

a2 


—_ Ls Vu + P) (32) 
dz 


The time-dependence represented by (28) is a result of 
the finiteness of the group velocity. It requires an inter- 
polation procedure, since the differentiations in the x and 
y directions require constancy of t as well as z, and unless 
the group velocity is constant, trajectories that are simul- 
taneous at one value of z will not be simultaneous at other 
values. Time-dependence may be of considerable importance in 
some applications, but there are many in which the scale of 
the problem is such that finite propagation velocity effects 
are not of interest. When this is the case, we may set 
v + @ in the above equations, and drop (28). Trajectories 
are then reinterpreted as space paths rather than space-time 
paths. If the material properties and boundary conditions 
are time-dependent, the trajectories will, of course, change 
with time, but no time derivatives will appear in the radia-~ 
tion propagation equations. 


A second simplification which is worth introducing is 
the restriction to axial symmetry. This is particularly con- 
venient in a developmental program, since it greatly reduces 
the number of trajectories required for a given spatial 


resolution. The time-independent axially symmetric forms of 
(31) and (32) may then be written 
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W, (2) = / 2mrI(r,z) dr (35) 
i 
is the beam power in the ring ry <r < Tia! and 


Titl 
a, (2) a 8 Qnra(r,z) I(r,z) dr/W, (2) (36) 


Tints 
ae 


is the effective absorption coefficient in the ring. 


It may be noted that (33) does not involve the radial 
variable r at all, and in (34) r appears as a dependent var- 
iable. This suggests that the right side of (34) be evaluated 
in terms of an independent variable which is related to W. 

It is also important to account properly for the radial 
boundary conditions, both at r = 0 and at the outer boundary. 
Fleck !®] devised a method of cylindrical cubic spline rep- 
resentation specifically for the r variable. An alternative 
approach becomes possible with other variable choices, which 
May provide certain advantages. 


For the evaluation of the right side of (34) we choose 
as dependent variable the cumulative area to radius r (per 


radian) 


and as 


where 


is the 
power. 


R-3564 
independent variable 
x = - log (1-S(r) /S,,) (38) 
r 
S(r) = f 2mrIi(r) dr (39) 
0 


beam power within radius r, and Si = S(~) is the total 
It is clear that a Gaussian beam profile 


2s 
I(r) = —S exp (mde Ar?) (40) 
TW 


is represented in these variables as the straight line 


y= ; wx (41) 


Since in the absence of perturbations the Gaussian 


profile propagates without change of shape (except for changes 


in the 


scale radius w), it is natural to assume a Gaussian 


Shape for the beam profile at the outer edge of the calcula- 


tional 


mesh as a boundary condition which is not likely to 


propagate spurious information into the interior. Such a 


boundary condition, in view of (41), is easily imposed with 


the variable definitions chosen. 


The derivative terms appearing in (26) and on the 


right side of (34) are easily evaluated. With derivatives 
by x denoted by dots, one finds from (37), (38) and (39) 


I = By a */Y (42) 
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ag ae. ae eee a ; | 

1) * 5F @ ~*~ Tay" (; : =) (43) | 
y y | 

d ee ae ee | 

TS, = dy Ty = ty = 2 et ey fyiey (44) 


d roe ROS ADE a Bahn ay eG «3 0¢6 
9" ay Tg = > FY + OFF = Se oer” as) 


Then (26) may be written 


4k Ir dr dr 
° 
1 € da 2 
= —> (= = (2yI T,) - yt 
4k I dy 1 L 
° 
ee 2 
== (27, + yT) + 2yT,) (46) 
4k 
° 
so that 
Che 3 2 
ee (47, + 07 © ont, + ayn.) . (47) 
° 


The refraction term in (34) is simply 


ou ” ee 
- ru/y (48) 


It may be noted that the right sides of (47), (48) and 
therefore of (34) vanish properly at r = 0. 


The derivatives of y and u with respect to x may be 
evaluated from spline representations. In order to construct 
these, trajectory values of x, y are first set up. These 


are determined by (37) and (38) following each integration 
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step for (33) and (34): 


where 


Se 


According to (45) derivatives of y up to fourth order are re- 
quired. A quintic spline representation provides continuous 


Cerivatives 


natural quintic spline code, QUINAT 


- log (1 - $,/8,), j= 1, n (49) 


0 (50) 


j 
= ow, Gata (51) 


i=l 
Set (52) 
2d 


through fourth order. A recently developed 


(9] was therefore modi- 


fied for the purpose. The procedure consists of two steps: 


es 


After an integration step for the Wy: Equation 
(33), and calculation of new values of a Equa- 
tion (49), the first part of the spline calcula- 
tion is done. This consists of setting up a 
matrix of coefficients and decomposing it into 
triangular and diagonal factors, using a 
Cholesky method. The results of this step de- 
pend only on the x. and can be used repeatedly 
for calculation of spline representation of 
other trajectory variables in step (2). 


The second step is done for each new or modified 
set of trajectory values of a, u and y. It 
yields coefficients Bye Ci, Di E , and Fs for 

J 


g 
the natural spline polynomials 


27 


R-3564 


3 


4 5 
+ E. (x-x.)° + F. (x-x. x, <x < x. 
j j pe By j+1’ 


j=l1,n (54) 


where the A are the prescribed values of some functions at 
the points e 5: The representation A(x) and its first four 
derivatives are continuous at the interior nodes x., j = 2, 
n- 1. At the end-points xy and Xn the third, fourth and 
fifth derivatives of A vanish automatically. 


(9] it was not possible 
to separate these two steps. By explicitly factorizing the 
coefficient matrix in step (1) this separation is made pos- 


In the original QUINAT program 


sible. Since the operations in step (2) must be carried out 
Many times compared to the more time-consuming ones in 


step (1), the gain in efficiency is large. 


It remains to specify integration methods for (33) 
and (34). The coefficients a. defined by (33) are evaluated 
from the spline representations of a(x) and y(x). With I(x) 
defined by (42) the integration in (36) can be done analyti- 
cally. If the a; vary slowly with z, they may be taken as 
constant within each increment of the z variable, in which 
case (33) may be integrated directly: 


W, (2 + Az) = W; (z) exp (- _ Az) (55) 


For 4 strongly dependent on z it may be preferable to inte- 
grate the equivalent form 


2. log W, = - a(z) (56) 


+ 


using a higher-order numerical procedure. 
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For integration of Equation (34) a predictor-correc- 
tor method is used. The method chosen is a generalization 
to nonuniform meshes of Stormer's method for solving second- 
order equations of the form 


5" Gie,z) . (57) 
dz 


Specifically, let 


h=2z - Zz 


The corrector, or closed, formula is then 


ey Pe 1 2 oat re 
Sr 41 = (h + g) ee h rn-1 + iz h (g” + gh h”) ony 
1 2 2 ‘e al 2 
2 = 1 2. ee 2 2, _v 
-g) ratl + 360 hg (g h') (2g + 5Sgh + 2h°)r (Zz, + ¢€) 


where the last term indicates the order of error incurred by 
dropping it. The predictor, or open formula chosen is 


% . 1 ee 
SF 44g ~ (et oie, core 6 nots a 
1 aoe 4 3.2 
ee h (h + g) (h + 29) rh * 37 (3g° h - 2h'g 
# hg) oo* te, + 2) (59) 


n 


Fo.iowing an initial estimate of ratl by (59), the right side 
of (34) is evaluated as described above to provide an esti- 

" 
mate of rel’ 
using the corrector (58). 


The process is then iterated to convergence | 
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The initialization of the trajectories does not pre- 
sent any particular problem. It is useful to record some 
relations characterizing Gaussian beams, in which the in- 
tensity varies as in (40) or (41). The scale radius w of 


such a beam has a z-dependence given by the hyperbola !10! 
2 2 4 2 
w(z)" =w+ (a= 2) (60) 
° kw" ° 
foe) 


where Zor w are the values of z and w at the waist of the 
beam. Each trajectory is a similar hyperbola 


Ge. (2. ) 
x; (2) = 1 w(z) (61) 
° 
so that 
a2 4r., 
=a r5 = = ri (62) 
dz kow (2) 


The result (62) is also deriverable from (34), (41) and (47) 
with the added assumptions a = p = 0. The values of r and 
r" needed to initialize trajectories in a Gaussian beam may 


be determined very simply by these relations. 


Material properties Kor a(r,z) and u(r,z) may be 
specified in a variety of ways. A dispersion relation for 
partially ionized gas, 


2 
<a, w 
Ko = 5a” oS oritusey ) (63) 
c & On 


where w_ is the electron plasma frequency and Vie V, are the 
electron collision frequencies with ions and neutrals res- 
pectively, is being used in test calculations. The 
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quantities a and u are obtained from (63) by extracting the 
complex square root and applying (19) and (20). 


A Fortran listing of the routines which perform the 
operations discussed above is presented in Appendix II. 


10. 
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APPENDIX I 


LISTING OF LSD WAVE PROPAGATION ROUTINES 
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